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ABSTRACT 

We present hydro dynamic simulations of the evolution of self-gravitating dense gas on scales of 
1 kiloparsec down to <parsec in a galactic disk, designed to study dense clump formation from giant 
molecular clouds (GMCs). These structures are expected to be the precursors to star clusters and this 
process may be the rate limiting step controling star formation rates in galactic systems as described 
by the Kennicutt-Schmidt relation. We follow the thermal evolution of the gas down to ^ 5 K using 
extinction-dependent heating and cooling functions. We do not yet include magnetic fields or localized 
stellar feedback, so the evolution of the GMCs and clumps is determined solely by self-gravity balanced 
by thermal and turbulent pressure support and the large scale galactic shear. While cloud structures 
and densities change significantly during the simulation, CMC virial parameters remain mostly above 
unity for time scales exceeding the free-fall time of GMCs indicating that energy from galactic shear 
and large-scale cloud motions continuously cascades down to and within the GMCs. We implement 
star formation at a slow, inefficient rate of 2% per local free-fall time, but even this yields global star 
formation rates that are about two orders of magnitude larger than the observed Kennicutt-Schmidt 
relation due to over-production of dense gas clumps. We expect a combination of magnetic support 
and localized stellar feedback is required to inhibit dense clump formation to ~1% of the rate that 
results from the nonmagnetic, zero- feedback limit. 

Subject headings: galaxies: ISM, galaxies: star clusters, methods: numerical, ISM: structure, ISM: 
clouds, stars: formation 



1. INTRODUCTION 

Star formation in galaxies involves a vast range 
of length and time scales, from the tens of kilo- 
parsec diameters and ~ 10^ yr orbits of galactic 
disks to the ~ 0.1 pc sizes and ~ 10^ yr dy- 
namical tim es of ind ividual prestellar cores (PSCs) 
(see ^cKee Ostrik er 2007; Scalo & Elmegreen 2004; 
[Elmegreen & Scalo 2004: iMac Low fc KlessenI l200l 
Ballesteros-Paredes et al., [20071 for reviews). Self- 
gravity in the gas is effectively countered by various 
forms of pressure support (including thermal, magnetic 
and turbulent), large-scale coherent motions (including 
galactic shear and large-scale turbulent ffows) that drive 
turbulent motions, and localized feedback from newborn 
stars in order to make the overall star formation rate 
relatively slow and inefficient at just a few percent con- 
version of gas to stars per lo cal dynamical timescal e 
across a wide range of densities ([Krumholz fc Tanll2007[ ). 
However, the relative importance of the above pro- 
cesses for suppressing star formation is unknown, even 
for the case of our own Galaxy. Other basic ques- 
tions such as "What is the typical lifetime of giant 
molecular clouds (GMCs)?", "What processes initiate 
star formation in localized clumps within GMCs?" , "Do 
star-forming clumps come close to achieving virial and 
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pressure balan ce?" and "Is the star cluster f ormation 
timescale lon^ (|Tan. Krumholz fc McKeell2006f ) or short 
fElme^greenl [20001 12QQ7D compared to free-fah?" are still 
debated. 

To investigate the star formation process within 
molecular clouds, a significant range of the internal 
structure of GMCs needs to be resolved including 
dense gas clumps expected to be the birth locations of 
star clusters. In one scenario of CMC formation and 
evolution, large-scale colliding atomic fiows have been 
invoked. Hi^h-resolution siniulations ([Folini fc W alderi 
T998[: iWalder fc Folinil [200(1 iKovama fc Inu tsuka 2001 
VazQuez-Semad eni et al.' '20031 '2011': 'van Loo etld] 
20071 . [2QlJ~ Hennebelle et al. 2008; Heitsch et al.] 
2008[: i Baneriee et al.i i2009l : lAudit fc Hennebellel 120101 : 
Clark et al. 2012, among others) show that cold, dense 
clumps and cores form in the swept-up gas of the shock 
front and that thermal and dynamical instabilities 
naturally give rise to a filamentary and turbulent cloud. 
However, there is little observational evidence that 
GMCs form from such large-scale and rapid converging 
fiows of atomic gas, especially in molecular-rich regions 
of galaxies, such as the Milky Way interior to the 
solar orbit. For example, even though GMCs are often 
associ ated with sp iral arms in galaxies, in the case of 
M51, [Koda et all (|2009[ ) find that the molecular to 
atomic gas mass fraction does not vary significantly 
from inter-arm to arm regions and infer relatively long 
CMC lifetimes ~ 100 Myr. The position angles of 
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projected angular momentum vectors of Galactic GMCs 
show ra ndom orientation s with respect to Galactic 
rotation (jKoda et all I2QQ6I : llmara Bl"it3l2Qll[ ). which 
is inconsistent wi th young, recently-f ormed GMCs in 
the simulations of iTasker fc T an (2009). Similar results 
are found even in relativ ely metal and m ol ecula r poor 
systems. In the LMC, iKawamura et al.l (|2009f ) infer 
GMC lifetimes of ~ 26 Myr, significantly longer than 
thei r dynamical t i mes. In M33. [RQSolowskv et al.l ()2003l ) 
and llmara et al.l ()2011f ) again find random orientations 
of the position angles of projected GMC rotation vectors. 

In this paper we investigate the processes of GMC evo- 
lution and star formation in a kiloparsec-scale patch of 
a galactic disk, starting from initial conditions in which 
GMCs have already formed. These are extracted from 
large, global galaxy simulations and then the evolution 
of the interstellar medium (ISM), especially GMCs, and 
star formation is followed over a relatively short timescale 
of ~ 10 Myr. This is less than the fiow crossing time 
of the simulation volume, but relatively long compared 
to the dynamical and free-fall times of the clouds. We 
achieve a minimum resolution of ~0.5 pc — enough to 
begin to resolve a significant range of the internal struc- 
ture of the GMCs. 

This paper, the first of a series, introduces the sim- 
ulation set up, and then investigates the effect of rel- 
atively simple physics, including: pure hydrodynamics 
(no magnetic fields, which are deferred to a future pa- 
per); a cooling function that approximates the transition 
from atomic to molecular gas and allows cooling all the 
way down to ~ 5 K; photoelectric heating; and simple 
recipes for star formation, parameterized to be a fixed 
formation efficiency per local free-fall time, e^. Local- 
ized feedback from newborn stars is difficult to resolve 
in this simulation set-up, so its treatment is deferred to 
a future paper. Thus the results of these simulations, 
which include the structural and kinematic properties of 
GMCs and clumps and their star formation rates, should 
be regarded as baseline calculations in a nonmagnetic, 
zero- feedback limit. As we shall see, by comparison with 
observed systems the degree of over production of dense 
gas and stars then informs us on the magnitude of sup- 
pression clump formation that is needed by these effects. 

2. METHODS AND NUMERICAL SET-UP 
2.1. Simulation Code and Initial Conditions 



The global galaxy simulations of iTasker fc TanI (|2009l 
hereafter TT09) followed the formation and evolution of 
thousands of GMCs in a Milky- Way- like disk with a ffat 
rotation curve. However, with a spatial resolution of 
~8 pc, only the general, global properties of the GMCs 
could be studied; not their internal structure. Following 
all of these GMCs to higher spatial resolution is very ex- 
pensive in terms of computational resources. Therefore 
we need to use an alternative method. 

Our approach is to extract a 1 kpc by 1 kpc patch 
of the disk, extending 1 kpc both above and below the 
midplane, centered at a radial distance of 4.25 kpc from 
the galactic center (see Fig. [1]). This is done at a time 
250 Myr after the beginning of the TT09 simulation, 
when the disk has fragmented into a relatively stable 
population of GMCs. We then follow the evolution of the 
ISM, especially the GMCs, including their interactions. 



internal dynamics and star formation activity, down to < 
parsec scales. These local simulations are able to reach 
higher densities, resolve smaller mass scales and include 
extra physics compared to the global simulations. 

Setting up the local simulations requires performing a 
velocity transformation to remove the circular velocity at 
each location and then adding a shearing velocity field 
so that the reference frame of the simulation is the local 
standard of rest at the center of the extracted patch of 
the galaxy disk. Periodic boundary conditions are intro- 
duced for the box faces perpendicular to the shear fiow 
(along the y-ax.is) and outfiow boundary conditions for 
the other faces. This set up is similar to shearing box 
simulations of astrophysical disks, but does not include 
the rotation of the frame of the patch (thus neglecting 
Coriolis forces). Since we only consider the evolution of 
the ISM over quite short timescales ~ 10 Myr (which is 
several local dynamical times of the GMCs, but much 
shorter than the orbital time of 130 Myr), this approxi- 
mation should not have a significant effect on the prop- 
erties of the clouds and their star formation activity. 

We also ignore radial gradients in the galactic po- 
tential, only including resulting forces perpendicular to 
the disk. Th e model for the potentia l is that used 
by TT09, i.e. iBinnev fc Tremaind (jl987[ ). evaluated at 
r = 4.25 kpc: 
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where Vc^o is the constant circular velocity in the limit 
of large radii, here set equal to 200 kms~^, Tc is the 
core radius set to 0.5 kpc, r and z are the radial and 
vertical coordinates respectively, and the axial ratio of 
the potential field is = 0.7. 

The grid resolution of the initial conditions is 128^ x 
256 which corresponds to a cell size of 7.8 pc and serves as 
the root grid for the high resolution simulations. Most of 
the simulations we present here involve 4 levels of adap- 
tive mesh refinement of the root grid, thus increasing the 
effective resolution to 2048^ x 4096 or about 0.5 pc. Re- 
finement of a cell occurs when the Jeans length drops 
below four c ell widths, in accorda nce with the criteria 
suggested bv lTruelove et all (|1997[ ) for resolving gravita- 
tional instabilities. 

The simulations performed in this paper were run using 
Enzo (jBrvan fc Normanlll997l : iBrvanll 19991 : 10 'Shea et all 
[200l . We use the second order Godunov scheme with 
the Local Lax-Friedrichs (LLF) solver and a piecewise- 
linear reconstruction to evolve the gas equations. Be- 
cause gas temperatures calculated from the total energy 
can become negative when the total energy is dominated 
by the kinetic energy, we also solve the non-conservative 
internal energy equation. We use the gas temperature 
value from the internal energy when the internal energy 
is less than one tenth of the total energy and from the to- 
tal energy otherwise. While this is higher than the often 
adopted ratio of 10~^, it does not affect the dynamics. 

2.2. Star Formation 

To model star formation, we allow collisionless star 
cluster particles, i.e. a point mass representing a star 
cluster or sub-cluster of mass M*, to form in our sim- 
ulations. These star cluster particles are created when 





Fig. 1. — View of the gas mass surface density, E^, over the 20 kpc diameter galactic disk of TT09 250 Myr after start of the simulation. 
The square is a 1 kpc sided region, enlarged in the right image, showing several GMCs. These are the initial conditions for the simulations 
of this paper. The intial clouds listed in Tableware marked. 



TABLE 1 
Star Particle Creation. 



Cell Size 
(pc) 


(cm~^) 


(yr) 


Min. Cell Mass 
(Mo) 


M*,niin 

(Mo) 


7.8 


100 


4.3 X 10^ 


1640 


1000" 


0.49 


10^ 


1.4 X 10^ 


400 


100 


0.125 


10^ 


4.3 X 10^ 


63 


10^ 



than 3000 K are prevented from forming stars. In fact, 
such conditions almost never arise in the simulations pre- 
sented here. 

When a cell reaches the threshold density, a star cluster 
particle is created whose mass is calculated by 



" Used in test runs not presented here and by Tasker (2011) 

^ Used in runs to be presented by Butler, Van Loo &; Tan, in prep. 

the density within a cell exceeds a fiducial star forma- 
tion threshold value of nH,sf = 10^ cm~^ for our 4-level 
refinement runs compared to a thres hold o f 10^ cm~^ in 
the lower resolution simulation of iTaskerl ()2011l ). (Note 
we assume nne = O.lnn so that the mass per H is 
2.34x 10~^^g). This density threshold is a free-parameter 
of our modeling, and its choice depends on the minimum 
mass that is allowed for star particles and the minimum 
cell resolution (see Table [T]). 

We use a prescription of a fixed star formation 
efficiency per local free-fall time, eff, for those re- 
gions with nn > ^H,sf- Relatively low and density- 
independent values of are implied by observational 
studies of GMCs (Zuckerman & Evans 197^) and their 
star- forming clumps (jKrumholz fc Tanll2007[ ). which mo- 
tivate our fiducial choice of 6ff = 0.02. Such values 
are also approximately consistent w ith numerical studies 
of tu rbulent, self-gravitating gas (jKrumholz fc McKed 
120051) and turbulent , self- gravitating, magnetized gas 
(jPadoan fc Nordhmdl [20Tll ) . Note that given the inabil- 
ity of our simulations to resolve individual star-forming 
cores, we do not impose requirements that the gas fiow 
be converging or that the gas structure be gravitationally 
bound in order for star formation to proceed. However, 
to rule out the possibility of star formation in the hot 
dense gas of shock fronts, cells with temperatures greater 



= eff^4 — At, 
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where p is the gas density, Ax^ the cell volume. At the 
numerical time step, and tff the free-fall time of gas in 
the ceh (evaluated as tff = (37r/32Gp)^/2). 

An additional computational requirement is the min- 
imum star cluster particle mass, M^^min, introduced to 
prevent the calculation from becoming prohibitively slow 
due to an extremely large number of low-mass parti- 
cles. If the calculated is < M^^min then a star clus- 
ter particle of mass M^^min is created with probability 
M*/M* ^min- In fact given the low value of e^, the fact 
that At <C tff and our adopted values of M^^min, we are 
always in this regime of "stochastic star formation" . 

The motions of the star cluster particles are calculated 
as a collisionless N-body system. Note these are not sink 
particles: there is no gain of mass by gas accretion. They 
interact gravitationally with the gas via a cloud-in-cell 
mapping of their positions onto the grid to produce a 
discretized density field. Note, however, that gravita- 
tional interactions between star cluster particles are thus 
softened to the resolution of the grid. In reality the dis- 
tribution of stellar mass represented by the star cluster 
particle would be spread out, but by amounts that are 
not set by the local gas density. Thus the structure of 
star clusters, made up of many simulation star cluster 
particles, is not well- modeled in our simulations and we 
do not present results on the details of the star clusters 
that form. 

2.3. Thermal processes 
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To describe the thermal behaviour of the ISM, we in- 
clude a net heating rate per unit volume given by 

H = nii[T - nuA] erg cm"^ s"\ (3) 

where F is the heating rate and A the cooling rate. Below 
we describe several different methods to model heating 
and cooling. 

The thermal processes introduce a new time scale, i.e. 
the cooling time tcooi = £^int/|^| where Eint =Pg/{l-^) 
is the internal energy, which is often much shorter than 
the dynamical time step set by the Courant condition. 
To prevent the cooling from increasing the numerical cost 
significantly, we set the numerical time step to the hy- 
drodynamical time step and sub-cycle the cooling with 
smaller time steps, i.e. O.ltcooi until the total cooling 
time step equals the dynamical time step. Then the tem- 
perature and internal energy are updated explicitly every 
sub-cycle after which the net heating rate and cooling 
time are recalculated. This sub-cychng prevents over- 
cooling and negative pressures by not resolving the cool- 
ing time. If the hydrodynamical time step is shorter than 
the cooling time, no sub-cycling is necessary. 

2.3.1. Simple Photoelectric Heating 

FUV radiation is absorbed by dust grains causing elec- 
trons to be discharged, which then heat the gas. This 
photoelectric (PE) heating has long been thought to be 
the dominant heating mechanism in the neutral ISM, in- 
cluding the relative ly low extinction portions of GMCs 
(|Wolfire et al .11 19951 ). We include a PE heating term 

FpE = 1.3 X IQ-^'^eGo erg s-\ (4) 

where e is the heating efficiency and G p the incident FUV 
field normalized to the lHabinj (|1968l ) estimate for the lo- 
cal ISM value. We set the heating effici ency to its maxi- 
muni value of e =0.05 for neutral grains ([Bakes fc TielensI 
[1994'). Also, we adopt Go = 4 as a value appropri- 
ate for the radial location of the simulate d region in a 
Milky Way- type galaxy, in agreement with IWoffire et al.l 
(|2003l ). Note that this estimated heating rate is approxi- 
mate in the sense that it does not follow the attenuation 
of the FUV field in the dense gas, nor its local generation 
from young star clusters. 

2.3.2. Simple Atomic Cooling 

We consider several cooling functions. Initially and for 
the higher temperature regime, we adop t the s olar metal- 
licity cooling curve of Sarazin & White (|1987l ) from tem- 
peratures of T = 10^ K down to T = 10^ K (although 
we do not expect gas temperatures in our simulations 
to exceed a few times 10^ K) and extend it d o wn to 
T = 300 K using rates from iRosen fc BregmanI (|1995l ) 
(see Fig. [2]). These temperatures take us to the upper 
end of the atomic cold neutral medium (|Woffire et al.l 
Il995l ). This cooling function is similar to that used 
by TT09. The floor temperature of 300 K acts as a 
minimum of thermal support against gravitational col- 
lapse and fragmentation. It imposes a minimum signal 
speed equal to the sound speed Cs = {"^kT / iirripY^'^ ^ 
1.80(T/300K)i/2 km s-\ where 7 = 5/3 and /i = 1.27 
(for an assumed nne = O.lnu). This signal speed is of 
the same order as observed velocity dispersions of clumps 
within GMCs fe.g. .Barnes et al...20111 
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Fig. 2. — Left: Cooling rate as function of temperature as de- 
fined by Sanchez-Salcedo et al. (2002) (solid) and Rosen & Breg- 
man with the high temperature range (> 10^ K) from Sarazin &; 
White (1987) (dashed). Right: The equilibrium temperatures for 
the Rosen & Bregman (dashed) and Sanchez-Salcedo et al. (solid) 
cooling curves assuming a constant PE heating term with Go = 4 
and for the Cloudy heating and cooling rates (squares). The dashed 
lines represent lines of constant temperature at 10^, 300 and 5 K. 
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Fig. 3. — Mean logarithmic visual extinction as a function of the 
density as derived using the Sanchez-Salcedo et al. simulation. The 
error bars show the dispersion on the mean, while the shaded area 
shows the distribution of all column extinctions in the numerical 
domain. The solid line represents the minimum column extinction 
due to absorption within the cell itself (see Eq. (6]). 
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However, the Rosen & Bregman cooling function does 
not include the formation and destruction of molecules or 
any cooling processes below the minimum temperature of 
300 K. By a combination of dust cooling and atomic and 
molecular line cooling, gas in real GMCs reaches temper- 
atures of ~ 5 — 10 K. We first take into account the effect 
of dust grains by adopting the atomic cooling function 
of lSanchez-Salcedo et alJ (120021), which mini ics the equi- 
librium phase curve of lWolfire et all ([19951 ) Figure [2] 
shows the difference between the Sanchez-Salcedo et al. 
and the Rosen & Bregman cooling rates. Note that, while 
the cooling rates are roughly the same for both curves, 
the Sanchez-Salcedo et al. cooling curve extends down 
to 5 K and has a thermally unstable temperature range 
between 313-6102 K. This gives rise to the co-existence 
of a cold and warm phase at the same pressure. 

2.3.3. Extinction- Dependent Heating and Cooling Functions 

A density- column extinction relation — The inclusion of 
molecular cooling is less straightforward since the for- 
mation of molecules depends strongly on the amount of 
attenuation of the radiation field. Molecules only form 
in the regions of relatively dense gas and dust that can 
shield their contents from destructive UV radiation. This 
means that the molecular cooling rate, and also heating 
rate, depends not only on density and temperature, but 
also on column extinction. 

Of these three variables, the column extinction is the 
only one that is not directly available during the simu- 
lation. For every grid cell an effective column extinction 
due to dust absorpt ion can be calculated usi ng a six-ray 
approximation (e.g. lGlover fc Mac-Lowl[2QQ7l ). Using the 
linear relation between column density and visual extinc- 
tion, i.e. Ay = 5.35 x 10~^^A/'h, the effective visual ex- 
tinction is given by 



2.5 



log 



1 ^ 

-^exp(-2.5Ay,i) 



mag. 



(5) 



Although the six-ray approximation already reduces the 
numerical cost of calculating the column extinction con- 
siderably, it remains time-consuming as it needs to be 
calculated for every cell at every time step. Still, an ad- 
ditional simplification can be made. In general, higher 
density gas has a higher column extinction than the sur- 
rounding lower density gas. We calculate the effective 
column extinction for the high-resolution run with the 
Sanchez-Salcedo et al. cooling function (see Sect. 13. 4p at 
10 Myr, but neglect cells within 10 grid spacings of the 
computational boundary. Figure [3] shows the full range 
of column extinctions in the numerical domain, i.e. the 
shaded area, and the mean logarithmic column extinc- 
tion as a function of the density with the error bars indi- 
cating the dispersion on the mean. We omitted density 
bins with fewer than ten cells contributing to the mean. 
While the extrema of the column extinction differ by 
more t han an order of m agnitude (similar to the results 
of e.g. iClark et"ani2012f ). the dispersion on the mean is 
much smaller. For example, at nn = 10 cm~^, 95% of all 
visual extinctions lie within the 0.458 ± 0.021 interval. 
Thus, the mean value is representative for the column 
extinction at a given density. An a posteriori check on 
the simulations using this relation (i.e. Run 5, 6 and 7 
of Table [2j) shows little deviation from the above result. 



The sudden increase at nn ^ lO^cm"^ is due to the 
resolution limitation, i.e. the effective visual extinction 
is dominated by absorption within the cell itself, i.e. 



Ax 

^y,eff = 5.35 X 10~^^^-nH mag 



(6) 



or Ay^eff = 4.1 mag for nu = 10^ cm~^ and Ax = 0.5 pc. 
The high density end is thus resolution-dependent and 
the column extinction is most likely overestimated. How- 
ever, for Ay > 10, the heating rate (and ionization bal- 
ance) is dominated by cosmic rays (see later), so this 
limitation should not affect the cooling and heating rates 
significantly. 

A table of heating and cooling rates — Using the density- 
column extinction relation we now generate a table of 
cooling and heating rates as function of density and 
temperature using the photodissociation code Cloudy 
(|Ferland et al.l[T998l version cOS.Ol). The routine is sim- 
ilar to the one of iSmith et al.l (|2008f ) and goes as follows. 
Firstly, for a given density, ranging from 10-^ to 10^ 
cm~^ with steps of 1 dex, the unextinguished local in- 
terstellar radiation field (jBlack 1987) with Go = 4 is 
directed through an absorbing slab. This slab has a con- 
stant density of nu = lcm~^, similar to the mean density 
of the disk in the initial conditions and a thickness corre- 
sponding to the visual extinction at the given density (see 
Fig. [3]). We assume the gas has solar abundances and in- 
clude a dust grain population with PAHs. The dust grain 
physic s use d in the code is de s cribed in Ivan Hoof et al.l 
(|2004f ) and IWeingartner et al.l (|2006f ). Background cos- 
mic rays are also included with a primary ionization rate 
of 2.5x10"^^ s~^, as well as the cosmic background ra- 
diation. 

The transmitted continuum, i.e. the sum of the at- 
tenuated incident and diffuse continua and lines, is then 
used as the radiation field incident on a gas parcel with 
the given density. The resulting cooling and heating 
rates, for a range of temperatures between 5 and 10^ K, 
are calculated self-consistently by simultaneously solv- 
ing for statistical and ionization equilibrium including 
all the necessary microphysics such as, among others, H2 
and CO formation and destruction. For temperatures 
above 10^ K we opt not to use Cloudy and simply adopt 
the cooling curve of iSarazin fc Whita ([1987) and set the 
heating rate to zero. From the generated table, the heat- 
ing and cooling rates for any density and temperature 
are derived using a bilinear interpolation. For densities 
and temperatures above and below the table limits, the 
rate of the limiting value is used. We implemented such 
a heating and cooling routine into Enzo. 

Figure [4] shows the calculated heating and cooling rates 
as a function of temperature for selected densities. The 
Sanchez-Salcedo et al. cooling curve is also plotted 
showing a close similarity in shape and magnitude for 
densities up to ~ 10^ cm~^. The deviation from the 
Sanchez-Salcedo et al. curve at higher densities stems 
from the onset of molecular cooling. Molecules form in 
regions with visual extinctions higher than Ay = 2.4 
or column extinctions highe r than 4.3 xlO^^ cm~^ (e.g. 
iTielens fc HollenbachI [19851 ) . Using the derived density- 
column extinction relation this translates to densities 
above ^ 10^ cm~^. Figure [5] shows the decomposition 
of the heating and cooling rates at thermal equilibrium 
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Fig. 4. — The cooling (left) and heating (right) rates as function 
of temperature for densities of lO*-* (short-dashed), lO"^ (dotted), 
10^ cm~^ (long-dashed) and 10^ cm~^ (dash-dotted) calculated 
using Cloudy. The cooling rates of Sanchez- Salcedo et al. (solid) 
are also plotted as reference. 
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log(nH) [cm~^] 

Fig. 5. — Cooling (solid) and heating (dashed) rates per unit 
volume as function of density for the Cloudy cooling function equi- 
librium temperatures given in Fig. [2] Only the processes that con- 
tribute the most are shown. 



(see Fig. [2] for the equilibrium temperatures). From this 
figure, we indeed find that molecular species, i.e. H2 and 
CO, contribute significantly to the heating and cooling 
above 10^ c m~^. Note that this figure is very similar 
to Fig. 8 of lGlover fc ClarkI (j2Q12[ ) in both the decom- 
position and the level of cooling rates. Note also that. 



the dust cooling disappears between ^ 10^ — 10^ cm~^. 
Up to 10^ cm~^, the dust grains are hotter than the 
gas. In these conditions, ion-grain collisions heat the 
gas due to thermal evaporation of neutralized ions (see 
Eq. 32 of [Baldwin et al.l (|1991[ )). while electron-grain 
collisions cool the gas. At densities below 10^ cm~^ the 
electron-grain collisions dominate the ion-grain collisions 
so that the net effect is gas cooling. As the density (and 
extinction) increases, less ionizing radiation can pene- 
trate the gas. Then ion-grain collisions start to domi- 
nate the electron-grain collisions and, consequently, the 
gas is heated. Above 10^ cm~^ the dust temperature falls 
below the gas temperature and grain collisions result in 
cooling of the gas. 

While the Cloudy cooling rates are adequately de- 
scribed by the Sanchez-Salcedo cooling function (up to 
10^ cm~^), the Cloudy heating rates deviate significantly 
from the PE heating rate we use with the atomic cooling 
functions, i.e. ^ 2.6 x 10~^^ erg s~^. It is clear that 
we have overestimated the heating efficiency by approxi- 
mately an order of magnitude. At high densities, the PE 
heating is no longer the dominant heating process (see 
Fig. [5j). Higher density regions tend to have higher dust 
extinctions of the external radiation field, thus blocking 
FUV penetration and PE heating. Then only cosmic rays 
can penetrate and heat the gas. The overall heating rate 
is further reduced by an order of magnitude. 

Not only does the decomposition of the cooling and 
heating rates give useful insights into the dominant pro- 
cesses at different temperatures and densities, it also pro- 
vides the opportunity to compare the numerical simu- 
lations with observations. For different emission lines, 
such as the 158 /im CII line and the 63 /im 01, tables 
similar to the cooling and heating rate tables can be 
constructed. Then emissivity maps and line profiles of 
optically-thin emission lines can be generated from the 
simulations during post-processing. Note, however, that 
such emission maps are only first order approximations 
as they do not include any radiative transfer, nor do they 
take into ac count the effects of time-dependent chemistry 
(e.g. iGlove r & Mac-Low 2007, 2011) or local abundance 
variations (^Shettv et al...2011a^ b). 



Thermal instability ?- 
curve (for nn < 10^ 



— Although the Cloudy cooling 
cm~^) has a similar temperature 
dependence between 300-6000 K as the Sanchez-Salcedo 
et al. curve, the equilibrium curve only exhibits a weak 
thermal instability (see Fig. [2j). Furthermore, the insta- 
bility is at lower densities than the Sanchez-Salcedo et al. 
curve. The shift of the instability towards lower densities 
is due to a higher c olumn extinctio n compared to the 10^^ 
cm~^ used by Wol fire et al.l (|1995f ). At low densities, the 
thermal equilibrium is set by the Lya cooling and PE 
heating. The increased attenuation reduces the electron 
fraction in the gas resulting in a lower PE heating (and 
thus a lower equilibrium temperature). 

Once the equilibrium temperature drops below 10^ 
K (and this happens at lower densities) , CII cooling 
becomes the dominant cooling process. I Wolfire et al.l 
(1995) show that the thermal instability can be at- 
tributed to CII cooling which causes a rapid drop in equi- 
librium temperature. The magnitude and presence of the 
thermal instability then depends on the temperature and 
density dependence of the CII cooling. Interstellar gas is 
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TABLE 2 
Set of simulations. 



Run 


AMR" 


Heating 


Cooling 


fi 


Star formation 


1 


No 


No 


RB^ 


1.27 


No 


2 


Yes 


No 


RB 


1.27 


No 


3 


Yes 


PE 


RB 


1.27 


No 


4 


Yes 


PE 


SS^ 


1.27 


No 


5 


Yes 


Cloudy 


Cloudy 


1.27 


No 


6 


Yes 


Cloudy 


Cloudy 


2.33 


No 


7 


Yes 


Cloudy 


Cloudy 


2.33 


Yes 



a- "No" implies min. resolution of 8 pc; "Yes" implies min. resolu- 
tion of 0.5 pc 

^ Rosen &; Bregman cooling function 

^ Sanchez-Salcedo et al. cooling function 

^ Photoelectric heating 

thermally unstable if (lFieldlll965[ ) 

nufdH\ fdH\ ^ 

By expressing the heating and cooling rate locally as 
a power- law of density and temperature, the instability 
constraint reduces to 

(a-6)-(a-/3)<l, (8) 

where a, resp. 6, is the density, resp. temperature, 
power-law index for the heating rate and a and j3 the 
indices for the cooling rate. We can use the above con- 
straint to understand the weak instability seen in Fig. O 
For temperatures between 300 — 6000 K the cooling rate 
has a power-law index of the order of 0.5, while the heat- 
ing rate shows no temperature dependence, i.e. {3 ^ 0.5 
and 6^0 (see Fig. [4j). Similarly, while the heating rate 
index, a, is roughly zero between 0.1 and 1 cm~^ (the 
range of densities for which the equilibrium temperatures 
is between 300 and 6000 K), the cooling rate decreases 
roughly as n^'^ . As a result, {a — h) — {a — j3) ^ 1. 
The instability criterion is thus only marginally satisfied 
explaining the weak thermal instability. 

3. RESULTS 

To understand the effect of different physical processes, 
we gradually include them in our models. All the details 
(e.g. included physics) of the simulations that we ran are 
listed in Table [1 

3.1. Measures of ISM Structure 

Several methods can be used to compare the different 
simulations and quantify the effect of the included phys- 
ical processes. We focus on measures that describe the 
changes in density and temperature. 

A visual analysis of the simulations is the simplest and 
quickest method of studying the geometrical changes. 
Column density plots show how the density structures 
change over time, with the maximum/minimum values 
indicative of increases/decreases in the density. 

We also follow the evolution of the mass fraction of the 
gas that is in defined density ranges: "GMC" gas is de- 
fined by nn > 10^ cm~^ and "Clump" gas is defined by 
riH > 10^ cm~^. Using our density-column extinction re- 
lation, this corresponds to a mean Ay^eff > 0.7 for CMC 
gas and Ay^eff > 40 for clump gas. While a visual extinc- 
tion of 0.7 does not imply molecular gas, it is only above 



TABLE 3 
List of the initial GMCs. 



Cloud 


Mass center 


Mass 


(7^ 








position (x,y,z) 


(106 Mo) 


(km s-l) 


(pc) 




A 


0.921, 0.164, -0.004 


0.09 


3.6 


17.9 


1.70 


B 


0.865, 0.249, 0.004 


0.06 


2.5 


15.1 


1.50 


C 


0.745, 0.232, 0.013 


7.47 


15.4 


47.9 


0.61 


D 


0.605, 0.597, 0.005 


2.53 


16.1 


37.4 


1.54 


E 


0.180, 0.159, -0.003 


0.79 


10.5 


21.5 


1.26 


F 


0.089, 0.588, -0.003 


1.00 


10.8 


27.0 


1.32 



" cr is the mass-weighted three-dimensional velocity dispersion. 
^ The radius is calculated from the cloud's volume assuming a 
spherical geometry. 

^ The virial parameter is calculated as oivir = ^g'm ' ^^^^re (7c is 
the one-dimensional velocity dispersion inside the cloud and given 
by cr2 = ^ + c2 (jPib et sl\\2Wl^ . 

nH = lO^cm"^ that we find gas with Ay^eff > 2.4 (see 
Fig. [3]). The density and column density structure can 
also be studied in more detail by using mass weighted 
probability density functions at different times. 

As the ISM has different phases, we also examine the 
mass fractions of the different temperature regimes, i.e. 
the cold, unstable, warm and hot phases. For the cold 
phase, T < 350 K. The unstable phase lies within the 
range of 350-6000 K. The warm phase extends from 6000 
to 10^ K after which the hot phase starts. These values 
are only indicative and we can even argue whether we 
need to include a thermally unstable range. Technically 
the Rosen & Bregman cooling function does not have 
one and the thermal instability in the Cloudy cooling 
function is weak. 

Finally, w e also use a clump-finding routine 

(jSmith et al.l I2009L for details) to identify individ- 
ual molecular clouds and derive their properties such 
as their mass, velocity dispersion and size. For this 
purpos e, we adapted the clump-finding routine included 
in YT (jTurk et al.ll2011[ ). This routine can also be used 
to identify dense clumps within the molecular clouds. 
However, as the highest resolution in our simulations is 
only 0.5 pc, clumps with sizes of the order of 1 pc are 
not properly resolved. Therefore, we refrain from any 
detailed analysis of the clumps. 

3.2. Initial Conditions 

Figure [1] shows the column density integrated along the 
z-axis and we can identify four distinct density struc- 
tures. In fact, the selected region contains six clouds. 
The two small clouds (A & B) are about to interact with 
a larger one (C) and are therefore not recognized as in- 
dividual clouds in the column density plot. Table [3] lists 
all of the clouds. Their properties span a large range in 
sizes and masses. The smallest cloud only has a radius 
of 15 pc and a mass of 6xlO^M0, while the largest cloud 
is a hundred times more massive (7.47 x 10^ Mq) and 30 
times larger in volume. The total mass in the clouds is 
1.2xlO^M0 which is about 70% of the gas mass in the 
simulation box. The clouds have diameters smaller than 
^ 100 pc, and thus have at most ~ 12 cells {Ax « 7.8 pc) 
across each linear dimension. To resolve turbulence in 
self-gravitating gas, iFederrath et al.l (|2011f ) suggest that 
the turbulent length-scale needs to be resolved by at least 
30 grid cells. The internal structure of the clouds and 
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their internal turbulence is thus not resolved in the ini- 
tial conditions. Note, however, that the velocity disper- 
sion of the clouds is significantly larger than the mini- 
mum sound speed of 1.8 km s~^ and increases with the 
cloud radius. Actually, the initial velocity dispersion of 
the clouds is high enough to give some approximate bal- 
ance against self-gravity. The mean virial parameter of 
the clouds is 1.32 with standard deviation of 0.35, close 
to the value of 1.3 with sta ndard deviatio n of .76 for 
Galactic GMCs derived by 'McK ee fc TanI (l2003h from 



analysis of the results of Solomo n" et al.l (19871). Note 
thejf CO-selected clouds studied bv lRoman-Duval et al.l 
(|2010l ). which trace somewhat higher densities, have me- 
dian virial parameters of 0.46. The values listed here are 
also in agreement with the values from the simulations 
oflDobbs et al. (2011). 

The temperature distribution of the initial conditions 
shows that 99% of the mass has a temperature below 
350 K and is therefore in the cold phase. All this gas 
lies within the galactic disk. The remaining 1% of the 
mass is diffuse hot gas with temperatures above 10^ K 
surrounding the disk. TT09 were not studying the full 
temperature structure of the warm gas and so did not yet 
include a heating term in their simulations. Thus most 
of the gas cooled down to the minimum temperature of 
300 K, with hotter components created in shocks. The 
coohng time for low density gas with temperatures above 
10^ K (i.e. the gas outside the disk) is of the order of 10"^ 
Myr and, thus, remains hot. As the hot gas lies outside 
the disk, its volume fraction can be used to derive a mean 
thickness of the disk, i.e. 140 pc. 



3.3. Evolution of Structure and Effect of Resolution 

We first carryout the simulation, running for 10 Myr, 
with the physics and resolution identical to the global 
galaxy simulation of TT09. 

For the given resolution o f ^ 7.8 pc and t he mi nimum 
temperature of 300 K, the iTruelove et al.l (|1997[ ) crite- 
rion, i.e. the Jeans length, Aj = (ttc^/Gp)-^/^ where G 
is the gravitational constant, associated with a grid cell 
should be at least 4 times larger than the size of the cell. 
Ax, is satisfied for densities up to nu 80cm~^. Similar 
considerations for the higher resolution simulations that 
include cooling down to ~ 5 K and reach much higher 
densities indicate that the Truelove criterion is not sat- 
isfied for resolving clump formation within GMCs. Arti- 
ficial fragmentation is expected to happen above densi- 
ties of 5x10^, 1.5x10^ and 2x10^ cm~^ for the Cloudy, 
Sanchez-Salcedo et al. and Rosen & Bregman cooling 
functions, respectively. 

Figure [6] shows the gas surface density for Run 1 (left 
panel) after 10 Myr. The clouds have not changed dra- 
matically over this time scale, but do show signs of grav- 
itational contraction, interaction (e.g. clouds A & B 
merge and collide with cloud C) and fragmentation (e.g. 
cloud D). The result of these interactions can be seen in 
the mass weighted PDFs (Fig.[7|), i.e. the "CMC" gas is 
redistributed over a slightly larger density range. How- 
ever, the same PDF also shows that the mass fraction of 
gas in "GMCs" is roughly the same as initially. In fact, 
the fraction remains quite constant throughout the sim- 
ulation, which spans a few free-fall times (Fig. [8]). The 
volume fraction of "CMC" gas shows an initial decrease. 



X (kpc) ™ 

Fig. 6. — Gas mass surface density along the 2;-axis after 10 Myr 
for the uniform Run 1 (left) and the AMR Run 2 (right). 
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Fig. 7. — Mass- weighted Probability Density Function (PDF) for 
Run 1 (red, dashed) and Run 2 (blue, dotted) after 10 Myr. We 
also plot the initial PDF (solid line). 

but also remains relatively constant. Thus the clouds 
are roughly in virial equilibrium as indicated by their 
initial and final virial parameters (see Tables [3] and [4j) . 
Thermal pressure and non-thermal motions within the 
clouds thus provide enough support against self-gravity 
to prevent runaway gravitational collapse (even though 
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Time [Myr] 

Fig. 8. — The mass fraction of gas in "GMCs" (nn > lO^cm"^; 
dashed) and in "Clumps" {nu > lO^cm"^; dotted) and the volume 
fraction of gas in "GMCs" (solid) for Run 1 (red) and Run 2 (blue). 



TABLE 4 

List of GMCs after 10 Myr for Run 1 and 2. 



Cloud 


Mass 


a 


R 




Run 1 


(106 Mo) 


(km s-i) 


(pc) 




A" 




















c 


7.7 


17.7 


44.7 


1.04 


Dl^ 


0.41 


5.1 


22.3 


1.10 


D2^ 


0.51 


4.8 


25.5 


0.73 


D3^ 


0.87 


6.9 


25.3 


0.92 


E 


0.70 


5.5 


26.1 


0.81 


F 


0.76 


6.7 


27.6 


0.99 




Cloud 


Mass 


a 


R 


Qlvir 


Run 2 


(106 Mo) 


(km s-i) 


(pc) 




A" 




















Cl^ 


1.1 


13.6 


21.9 


2.2 


C2^ 


6.9 


25.6 


29.0 


4.1 


Dl^ 


0.10 


6.4 


9.3 


3.1 


D2^ 


0.12 


5.5 


10.2 


2.0 


D3^ 


0.48 


16.4 


12.3 


5.4 


D4^ 


0.12 


14.8 


18.0 


1.9 


D5^ 


0.19 


13.0 


6.4 


4.5 


D6^ 


0.11 


8.4 


6.4 


2.9 


E 


0.70 


15.2 


14.7 


2.9 


F1& 


0.16 


6.8 


9.6 


1.6 


P2b 


0.81 


12.4 


15.4 


2.0 



" Cloud A and B merge with cloud C. 

^ The cloud fragments into multiple clouds. 



the non-thermal motions are not well resolved in this uni- 
form grid simulation). A similar observation was made 
by TT09 in their global simulation, i.e. the gas properties 
in the full disk are quasi-steady after timescales of about 
150 Myr. However, note that now the velocity frame of 
the simulation volume is the local standard of rest of the 
disk at the center of the box, so the fast, ~ 200km s~^ or- 
bital velocities that were present in the TT09 simulation 
are now absent. Modeling these fast circular velocities on 
a finite rectilinear grid led to relatively large numerical 
viscous heating that helped stabilize the TT09 GMCs. 
Thus it is not surprising that we now see that the clouds 
are able to evolve to somewhat higher densities than were 
seen by TT09. 



By increasing the resolution up to ~ 0.5 pc the GMCs 
can now be better resolved and the evolution of dense 
clumps within the clouds begins to be captured. While 
the "CMC" mass fraction increases slightly to 72%, the 
volume fraction decreases by a factor of 2 within 2 Myr 
after which it remains roughly constant (see Fig. [8]). In 
approximately the same time span, about half of the 
"CMC" gas accumulates in "Clumps" (see Fig. [8]). Then 
the mass fraction in clumps remains nearly constant and 
reaches a new quasi- steady state. So, initially, the gas 
within the clouds collapses to form filaments with dense 
clumps due to the increased resolution. The clouds thus 
contract, but, at the same time, the velocity dispersion in 
the clouds increases. Thermal pressure and non-thermal 
motions again counter the effects of self-gravity to viri- 
alize the cloud (the virial parameters of the clouds are 
given in Table [4]). The initial evolution of the clouds is 
thus a direct consequence of the increase in resolution. 

The formation of clumps can also be seen in the PDFs 
(see Fig.[7j). While the PDFs for Run 1 and 2 are similar 
up to Tin = 80 cm~^, the gas above the "CMC" density 
threshold is redistributed towards higher densities. The 
higher densities also give rise to higher surface densities, 
i.e. the maximum value increases to ~ 10^ g cm~^ (see 
Fig.©. 

3.4. A Multiphase Interstellar Medium 

While the increased resolution helps to describe the 
substructures of GMCs in greater detail, the thermal 
properties of the ISM are poorly reproduced. As only 
cooling is included in Run 1 and 2, most of the gas 
within the disk is at the fioor temperature of 300 K. 
To reproduce the multiphase character of the ISM, i.e. a 
cold, dense and a warm, diffuse phase, we include diffuse 
heating. Additionally, we study the infiuence of differ- 
ent cooling functions. While atomic cooling can be ade- 
quately described by the Rosen & Bregman cooling func- 
tion (used in Sect. 13. 3|) . the Sanchez-Salcedo et al. func- 
tion extends down to a temperature of 5 K, i.e. nearly 
two orders of magnitude lower, and includes a thermally 
unstable temperature range. Extinction-dependent cool- 
ing, especially from CO molecules, is only taken into ac- 
count in the Cloudy cooling function. 

Including diffuse heating mostly affects the gas outside 
the GMCs. Figure [2] shows that, for nn < 90 cm~^, the 
gas has a higher equilibrium temperature than 300 K. For 
the other cooling functions the critical density for heating 
is at lower densities, i.e. 10 cm~^ for Sanchez-Salcedo et 
al. and 1 cm~^ for Cloudy. The diffuse intercloud gas 
thus moves from the cold phase to the warm phase. The 
decrease in the cold gas mass fraction is most significant 
for the Rosen & Bregman cooling function. Only the 
"CMC" gas which is 70% of the total, remains in the 
cold phase. For the Sanchez-Salcedo et al., resp. Cloudy, 
cooling function, some of the gas surrounding the clouds 
is also in the cold phase so that the mass fractions are 
85%, resp. 96%. The mass fraction of gas that is in the 
cold phase (i.e. molecular gas and CNM) in the inner 
Galaxy disk region of the molecular ring is ~ 0.5 (Wolfire 
et al. 2003). The simulation values for this mass fraction 
are somewhat higher, which we expect is due mostly to 
their present lack of ionization and supernova feedback 
processes. 

Together with the temperature, the pressure in the in- 
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Fig. 9. — Same as Fig. 7, but for Run 2 (green, dash-dotted), 
Run 3 (solid). Run 4 (red, dashed) and Run 5 (blue, dotted). 
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Fig. 10. — Same as Fig. 8, but for Run 3 (black). Run 4 (red) and 
Run 5 (blue). We also plot the volume fraction for Run 6 (green). 
The mass fraction for Run 6 is plotted in Fig. [15] 



ter cloud region increases significantly. For example, for 
nn = 1 cm~^, the pressure difference is of the order 
lO^/c^, where /cb is the Boltzmann constant. The asso- 
ciated heating time scales are of the order of 0.1 Myr 
for 0.1 cm~^, but are much shorter and even below the 
numerical time step for nn > 1 cm~^. (Remember that 
the time step is determined by the Courant condition 
for the hydrodynamics and not limited by the cooling 
time.) As a consequence, a significant amount of energy 
is added during the first time step to the simulation, i.e. 
of the order of 10^^ erg for the Rosen & Bregman and 
the Sanchez-Salcedo et al. functions and 10^^ erg for 
the Cloudy function. Note that this initial adjustment 
is unphysical, and can be regarded as a transient asso- 
ciated with the initial conditions. The added energy is 
primarily deposited near the midplane of the disk and 
eventually causes the disk to expand. The mean disk 
thickness for Run 3 and Run 4 (the atomic cooling func- 
tions) is ~ 600 pc, an increase of a factor of 4, after 10 
Myr. Because of the lower energy deposit for the Cloudy 
function, the mean disk thickness of Run 5 increases only 
by 30% to 190 pc. The larger disk of Run 3 and 4 can 
also be seen in Fig. [9] where the PDFs show that the den- 
sities between 10~^ and 1 cm~^ contain more mass than 
in the simulations without heating. 

While the higher external pressure causes the disk to 
infiate, it also acts as an additional force confining the 
molecular clouds. The higher external pressure resulting 
from this heating of the disk is, however, much smaller 
than the internal pressure of the GMCs, which is set by 
their self-gravitating weight. 

For the runs with the Sanchez-Salcedo et al. and 
Cloudy cooling functions, the gas in the interior of the 
clouds cools and loses pressure. This, potentially, has a 
large effect on the cloud evolution. However, by compar- 
ing the volume fraction of "CMC" gas in Run 4 and 5 
(when the cloud loses internal thermal pressure) with the 
one in Run 3 (where the internal thermal pressure of the 
cloud stays constant), we find that the resulting effect is 
minimal (see Fig. [TQ|) . The mass distributions above 10^ 
cm~^ are nearly identical with a small increase to 75% 
for Run 4 and to 79% for Run 5 (see Figs. [10] and [9]). 
Also, the amount of gas that ends up in dense clumps 
is independent of the cooling function. Furthermore, the 
clouds found in Run 3, 4 and 5 after 10 Myr using the 



cloud-finding algorithm have similar masses, velocity dis- 
persions and sizes. 

Since the gas in the clouds is cold and predominantly 
molecular, not atomic, we also ran a simulation with the 
Cloudy cooling function and where we use \i = 2.33 in- 
stead of /i = 1.27. This change does not affect the dy- 
namics of the clouds as can be seen in Figs. [TOl [13] and 

[m 

While the global properties of the clouds are simi- 
lar, the cloud substructure, i.e. the clump distribution, 
changes with the cooling function (see Fig. [TT]) . Much 
more filamentary and clumpy structures are present for 
the Cloudy cooling function than for the other two. How- 
ever, this is partly a numerical effect due to the applied 
refinement criterion and to the timestep used for evolving 
the simulation. As the different cooling functions tend 
to cool the gas to different equilibrium temperatures (see 
Fig. [2]), the densities at which a grid cell is refined differ 
as the gas temperature influences the Jeans length. Such 
changes introduce s mall variations in AMR simulations 
(|Niklaus et al.ll2009f ). 

3.5. Star Formation 

The high resolution simulations described above (Runs 
2-6) evolve to a state where a large fraction, 40-50%, of 
the gas is in clumps, as deflned by nu > 10^ cm~^. Now, 
in Run 7, we introduce star formation in these objects, 
following the method described in §2.2. As the density 
threshold for gas in dense clumps is the same as the crit- 
ical density of our star formation routine, star particles 
will form in the dense clumps. 

Figure [12] shows the gas surface density of Run 7 with 
the star cluster particles plotted on top. The star cluster 
particles are concentrated within the molecular clouds. 
Note, however, that star particles can be ejected from the 
clouds, especially in clouds with large amounts of angu- 
lar momentum, e.g. from a collision. The star formation 
has not changed much of the global density structure 
or dynamics, but has reduced the maximum gas surface 
densities by about an order of magnitude compared to 
Run 6 without star formation. This can also be seen in 
the PDFs where the star formation process only produces 
a deviation above nn ^ 10^ cm~^ and limits the maxi- 
mum density in the simulation to ~ 10^ cm~^ (Fig. [T3|). 
As we do not yet include stellar feedback in our simu- 
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Fig. 12. — Same as Fig. 6 but for Run 7. The white dots represent 
the star particles. 



Fig. 11. — Same as Fig. 6 but for Run 3 (top left), Run 4 (top 
right) and Run 5 (bottom). 



lations, the stars only interact gravitationally with their 
maternal cloud. 

Figure [TH shows a rendered visualization of the ISM 
structures at the end of Run 7, highlighting the volume 
density thresholds that define "GMCs" and "Clumps". 
Several hundred parsec long filaments of high density gas 
are visible, which bear a qualitative resemblance to some 
observed Galactic infrared dark cloud structures (e.g. the 
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Fig. 13. — Mass-weighted PDF for Run 6 (red, dashed) and Run 
7 (blue, dotted). We included the stellar mass in the total mass to 
normalize the distribution function. 




Fig. 14. — Volume rendered number density of Run 7. The 
"GMC" threshold volume density, nn ~ 100 cm~^, is coloured 
blue, while the "Clump" threshold volume density, nn ~ lO^cm"^, 
gas is coloured red. 



"Nessie" nebula. ^Jack son et al.ll2QlQf ). A comparison of 
dynamical st ate of the simulation fila ments with observed 
IRDCs (e.g. [Hernandez et al.ll2Q12f ). will be carried out 
in a subsequent paper. 

The removal of mass from the gas phase because of star 
formation can be seen by following the gas mass fraction 
in molecular clouds and in dense clumps for Run 6 and 
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Fig. 15. — Same as Fig. 8, but for Run 6 (red) and Run 7 (blue). 
The stellar mass fraction for Run 7 (solid) is also shown. For Run 7 
the normalization is done with the sum of the total gas and stellar 
mass. 
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Fig. 16. — Left: The star formation rate (Ssfr) as function 
of the gas surface density. The symbols show the observations 
of Bigiel et al. (2008), while the solid line shows the evolution of 
^SFR in Run 7. The arrow shows the direction of the evolution. 
The star formation rate quickly rises to more than 100 times the 
observed value — a result both of the initial condition having a 
relatively high mass fractions of dense gas and the lack of support 
from magnetic fields or disruption by stellar feedback. Gas is con- 
sumed, but the SFR remains about 100 times larger than the levels 
seen in galactic disks with similar gas content. Right: XlgPR (solid) 
as function of time for Run 7. The dashed line shows the rescaled 
mass fraction of gas in dense clumps for Run 7. 



7 (Fig. [15]). After 2 Myr, dense clumps start to form 
and immediately produce stars. As molecular gas is con- 
verted into stars, the mass fraction of gas in molecular 
clouds starts to decrease (Note that we calculate the mass 
fraction to the total mass including the stellar mass). Af- 
ter 10 Myr, the mass fraction in "GMCs" in Run 10 is 
about 60% lower than in Run 6. Similarly, the dense 
clump mass fraction drops to 2% of the total mass com- 
pared to 32% in Run 6. As the stellar mass is about half 
of the total mass after 10 Myr and star cluster particles 
only originate within dense clumps, this suggests that 
the dense clump gas is continuously replenished from the 
lower density molecular gas. This replenishment, how- 
ever, does not keep up with the rate at which the dense 
clumps convert gas into stars. The free-fall time of gas 
within the clump is shorter than the free-fall time of the 
region surrounding the clump, i.e. the mean density de- 
creases when including the gas around the clump. Hence, 
the mass fraction in dense clumps decreases with time. 
As the star formation rate depends on the gas mass in 
dense clumps (Eq. [2|), the evolution of the star forma- 
tion rate should follow the curve of the mass fraction of 
the dense clumps. Figure [16] indeed shows that this is 
the case with a maximum star formation rate of 1.8 Mq 
yr~^ kpc~^ after 2 Myr and then a steady decrease to 0.2 
M0 yr~^ kpc~^ at 10 Myr. This is more than 2 orders of 
mags:n itude more than the star formation rates observed 
by Bi giel et al.l (j2008[ ) (see Fig. [T6]) . 

Of course, the star formation rates in Run 7 depend 
on the values of the parameters and potentially on 
M^^min and nH,sf in the star formation routine. We did 
additional simulations where we varied these parameters. 
Increasing MH<,min to 200 Mq shortens the run time of the 
simulation as the number of star particles in the simu- 
lation decreases, but it does not change the star forma- 
tion rate. The total stellar mass only deviated by less 
than 2%. By changing nH,sf to 10^ cm~^, stars start 
to form earlier as these densities are reached at earlier 
times. However, the overall star formation rate is not 
affected very much: the total stellar mass only increases 
by 6%. On the other hand, if we increase nH,sf to 10^ 
cm~^ the stellar mass decreases by ^ 26%. The increased 
critical density reduces the gas reservoir from which the 
star particles form, although only by a relatively small 
amount. 

4. CONCLUSION AND DISCUSSION 

We have investigated the star formation process down 
to <parsec scales in a galactic disk. We extracted a 
kiloparsec-scale patch of the disk from the large, global 
simulation of TT09 and increased the resolution down 
to 0.5 pc. This allowed us to study the structure and 
evolution of GMCs in greater detail. We also included 
additional physics such as heating, atomic and molecular 
cooling, and a simplified approach to star formation. So 
far we have neglected the countering effects of magnetic 
fields and localized stellar feedback. 

We used a novel approach to include molecular cool- 
ing in our models. The formation of molecules de- 
pends strongly on the amount of attenuation of the ra- 
diation field. From a hi^h-resolu tion simulation includ- 
i ng th e atomic cooling function of lSanchez-Salcedo et al.l 
(I2OO2D that reprodu ces the equilibrium phase curve of 
iWolfire et al.l (|l995h we find that the column extinction 
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can be expressed as a function of gas density. Such a one- 
to-one relation eliminates the need for time-consuming 
column-extinction calculations to assess the attenuation 
of the radiation field in the numerical simulations. We 
also use this extinction-density relation to generate a ta- 
ble of cooling and heating rates as function of density and 
temperature with the code Cloudy. The resulting cool- 
ing function resembles the atomic cooling function up 
to densities of 10^ cm~^, above which molecular species 
start to dominate the cooling rates. However, we need to 
keep in mind that the heating and cooling rates are only 
first order approximations as the extinction law is only 
a mean relation and as local abundance variations and 
time-dependent chemistry are not considered. Further- 
more the simulations do not take into account the local 
generation of FUV radiation from young star clusters. 

With an increased resolution of ~ 0.5 pc our simula- 
tions are able to capture a significant range of the internal 
structure of molecular clouds. While the global proper- 
ties, such as the mass in the molecular clouds, remain the 
same, filaments and dense clumps form within the clouds 
shifting the mass distribution towards higher densities. 
The mass distribution within the molecular clouds are 
independent of the applied cooling function even though 
the three cooling functions describe different aspects of 
the thermal properties of the ISM. This suggests that 
the thermal pressure is of minor importance within the 
gravitationally-bound clouds. Then self-gravity and non- 
thermal motions determine the cloud structure. The 
nonthermal motions are driven by bulk cloud motions 
inherited by the GMCs that were formed and evolved in 
a shearing galactic disk where cloud-cloud col lisions are 
freque nt and infiuence GMC dynamics (TTQQ. pTan et al] 
(|2Q12f )). Our simulations begin to resolve the cascade of 
the kinetic energy from these larger scales and processes 
down to the smaller scales of clumps in the GMCs. This 
approach is to be contrasted with the method of driv- 
ing turbule nce artificially in per iodic box simulations of 
GMCs fe.g. lSchmidt et al.l[2QTQl ) or of forming GMCs in 
large-scale converging fiows of atomic gas where the tur- 
bule nce is driven by thermal and dynamical ins tabilities 
(e.g. iHeitsch et aIlf2QQ8l : IVazquez-Semadeni et al. 2011 ) . 

Of the gas within molecular clouds 50-60% is in dense 
clumps with nu > 10^ cm~^. This value is much higher 
than observed in nearby GMCs, e.g. 90% of the clouds 
in the Bolocam Galactic Plane Survey have a ratio of 
clump mass to cloud mass, or clump fo rmation efficiency, 
between and 0.15 (jEden et al.ll2012[ ). The high clump 
formation efficiency is partly due to our resolution limit. 
We do not properly capture the formation of individ- 
ual PSCs. so that the turbulent dissipation range is not 
fully resolved. The clumps then lack turbulent support 
against self-gravity hereby attaining higher densities and 
accumulating more mass. 

The surface densities of the clumps in Run 7 are in 



the ran^e o f ~0.1 t o ^10 g cm~^. Galactic IRDCs (e.g. 
Butler fc T an 2009, l2012f ) and star-forming clumps (e.g. 
Mueller et al.. ■2002f ) are found in the range ~0.1 to a 
few g cm~^. The most extreme mass surface density 
clumps seen in our simulations are probably prevented 
from forming in reality by localized feedback processes 
from star formation, espe cially the momentum input 
from protostellar outffows (jNakamura fc Lil [20071 ). 

The star formation rate in our simulations exceeds 
that expected from the Kennicutt- Schmidt relation (e.g. 
iBigiel et al.l l2011f ). These authors find a close rela- 
tion between the surface density of H2 and the star 
formation rate surface density for 1 kpc resolution re- 
gions within 31 disk galaxies. With a surface density 
of roughly 8 Mq pc~^ in our simulations, the star for- 
mation rate of 0.3 M© yr~^ kpc~^ after 10 Myr is 
about 100 times higher than observed. This over ef- 
ficiency of star formation in our simulations is a sim- 
ple refiection of the high mass fraction in dense clumps: 
while the GMCs are globally stable, there is no support 
against free-fall collapse in local regions of the GMCs. 
We can identify two physical mechanisms of reducing 
this n iass fraction: magnetic fi e lds an d stellar feedback. 
From iChandrasekhar fc Fermil ()1953l ) we can estimate 
the critical magnetic field required for support against 
self-gravity (neglecting the contribution of thermal and 
turbulent support), i.e. 

B.rit = 27rRpVG. (9) 

Assuming that the clumps in a GMC condense out of 
local volumes of radius ~ 10 pc and a GMC density of 
Tin = 100 cm~^, we find that a mean magnetic field of ~ 
10 /iG is sufficient. This is si milar to the observatio nally 
inferred value at this density (|Crutcher et al.|[2010[ ). We 
will study the effect of magnetic support in a subsequent 
paper. 
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